home *** CD-ROM | disk | FTP | other *** search
/ CD-ROM Today - The Disc! 5 / CD-ROM Today - The Disc (Issue 5)(November 1994).ISO / mac / Mac shareware / Education / RLaB / toolbox / lagrange.r < prev    next >
Encoding:
Text File  |  1994-09-21  |  1.6 KB  |  65 lines  |  [TEXT/ttxt]
  1. //-------------------------------------------------------------------//
  2. //LAGRANGE  Lagrange interpolation of arbitrary order. Y = LAGRANGE(TAB,X0,N)
  3. //          returns an N-th order interpolated value from table TAB, looking
  4. //          up X0 in the first column of TAB.
  5.  
  6. //          NOTE: TAB's 1st column is checked for monotonicity. It is an
  7. //          error to request a value outside the range of the first column
  8. //          of TAB for X0.
  9.  
  10. //          Original Author: Michael F. Saucier 10-16-87
  11. //-------------------------------------------------------------------//
  12.  
  13. lagrange = function (tab, x0, N)
  14. {
  15.  
  16.   local (dx, i, j, jj, jmin, lden, lnum, m, n, seq, sig, tab2, y0);
  17.  
  18.   if (!exist (N)) { 
  19.     error("Wrong number of input arguments."); 
  20.   }
  21.  
  22.   dx = diff (tab[;1]);
  23.   sig = sign (dx[1]);
  24.   if (any (sign (dx) - sig))
  25.   {
  26.     error("First column of the table must be monotonic.");
  27.   }
  28.  
  29.   // Check range of 1st column versus x0
  30.  
  31.   if (tab[1;1] > x0 || tab[tab.nr;1] < x0) {
  32.     error ("lagrange: x0 must be within range of tab");
  33.   }
  34.  
  35.   i = find (tab[;1] == x0);
  36.   if (i != 0) 
  37.   {
  38.     return tab[i;2];
  39.   }
  40.  
  41.   m = size (tab)[1]; n = size (tab)[2];
  42.   jmin = min ([max ([min (find (tab[;1] > x0)) - fix ((N+1)/2),1]), m - N]);
  43.   tab2 = tab[jmin:jmin+N;];
  44.   jj = 1:N+1;
  45.  
  46.   seq = x0*ones (1, N+1) - tab2[jj;1]';
  47.   lnum = prod (seq) ./ seq;
  48.  
  49.   lden = ones (1, N+1);
  50.   for (i in jj)
  51.   {
  52.     for (j in jj)
  53.     {
  54.       if (j != i) 
  55.       {
  56.     # fprintf ("stderr", "i = %i, j = %i\n", i, j);
  57.         lden[i] = lden[i] * (tab2[i;1] - tab2[j;1]);
  58.       }
  59.     }
  60.   }
  61.  
  62.   y0 = sum (lnum' ./ lden' .* tab2[jj;2]);
  63.   return y0;
  64. };
  65.